Processing Pipeline for Hematopoesis Data from Popescu et al (2019)

Bacher Group

Author
Affiliation

Jack R. Leary

University of Florida

Published

October 11, 2024

1 Libraries

Before we do anything else we load in the necessary software packages in both R & Python.

1.1 R

Code 
library(dplyr)       # data manipulation
library(Seurat)      # scRNAseq tools
library(scLANE)      # trajectory DE
library(ggplot2)     # pretty plots
library(biomaRt)     # gene annotation
library(tradeSeq)    # more trajectory DE
library(patchwork)   # plot alignment
library(slingshot)   # pseudotime estimation
library(reticulate)  # Python interface
select <- dplyr::select
rename <- dplyr::rename

1.2 Python

Code 
import scvi                                                     # integration
import warnings                                                 # filter out warnings
import numpy as np                                              # linear algebra tools
import scanpy as sc                                             # scRNA-seq tools
import pandas as pd                                             # DataFrames
import anndata as ad                                            # scRNA-seq data structures
import scvelo as scv                                            # RNA velocity
import cellrank as cr                                           # cell fate estimation
import matplotlib.pyplot as plt                                 # pretty plots
from scipy.io import mmread, mmwrite                            # sparse matrix IO
from cellrank.estimators import GPCCA                           # CellRank estimator
from scipy.sparse import coo_matrix, csr_matrix                 # sparse matrices
from cellrank.kernels import PseudotimeKernel, CytoTRACEKernel  # CellRank kernels

Here we set a few global variables to reduce the amount of printed output caused by our Python code.

Code 
warnings.simplefilter('ignore', category=UserWarning)
sc.settings.verbosity = 0
scv.settings.verbosity = 0
cr.settings.verbosity = 0

2 Visualization tools

To make our visualizations prettier we’ll define some consistent color palettes, and set a theme for matplotlib that matches ggplot2::theme_classic().

2.1 Color palettes

Code 
palette_heatmap <- paletteer::paletteer_d("MetBrewer::Cassatt1", direction = -1)
palette_celltype <- as.character(paletteer::paletteer_d("ggsci::category10_d3"))[1:6]
names(palette_celltype) <- c("HSC", "Kupffer cell", "Monocyte-macrophage", "Monocyte", "Monocyte precursor", "Neutrophil-myeloid progenitor")
palette_cluster <- paletteer::paletteer_d("ggsci::default_igv")

2.2 Theme for matplotlib

Code 
base_size = 12
plt.rcParams.update({
    # font
    'font.size': base_size, 
    'font.weight': 'normal',
    # figure
    'figure.dpi': 320, 
    'figure.edgecolor': 'white', 
    'figure.facecolor': 'white', 
    'figure.figsize': (6, 4), 
    'figure.constrained_layout.use': True,
    # axes
    'axes.edgecolor': 'black',
    'axes.grid': False,
    'axes.labelpad': 2.75,
    'axes.labelsize': base_size * 0.8,
    'axes.linewidth': 1.5,
    'axes.spines.right': False,
    'axes.spines.top': False,
    'axes.titlelocation': 'left',
    'axes.titlepad': 11,
    'axes.titlesize': base_size,
    'axes.titleweight': 'normal',
    'axes.xmargin': 0.1, 
    'axes.ymargin': 0.1, 
    # legend
    'legend.borderaxespad': 1,
    'legend.borderpad': 0.5,
    'legend.columnspacing': 2,
    'legend.fontsize': base_size * 0.8,
    'legend.frameon': False,
    'legend.handleheight': 1,
    'legend.handlelength': 1.2,
    'legend.labelspacing': 1,
    'legend.title_fontsize': base_size, 
    'legend.markerscale': 1.5
})

3 Helper functions

We first define a utility function to make our plot legends cleaner to read & easier to make.

Code 
guide_umap <- function(key.size = 4) {
  ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(size = key.size,
                                                                    alpha = 1, 
                                                                    stroke = 0.25)))
}

4 Data

4.1 Import hematopoesis data

We start by importing the AnnData object we downloaded from the Human Developmental Cell Atlas portal; it contains a variety of celltypes beyond just the myeloid lineage which we’re interested.

Code 
ad_liver = ad.read_h5ad('../../data/hematopoesis/fetal_liver_alladata.h5ad')

We subset to just the myeloid lineage, and remove cells with a percentage mitochondrial reads greater than 15%.

Code 
ad_myeloid = ad_liver[ad_liver.obs['cell.labels'].isin(['HSC_MPP', 'Neutrophil-myeloid progenitor', 'Monocyte precursor', 'Monocyte', 'Mono-Mac', 'Kupffer Cell'])]
ad_myeloid.obs['celltype'] = (
    ad_myeloid.obs['cell.labels']
              .map(lambda x: {'HSC_MPP': 'HSC', 'Mono-Mac': 'Monocyte-macrophage', 'Kupffer Cell': 'Kupffer cell'}.get(x, x))
              .astype('category')
)
ad_myeloid = ad_myeloid[ad_myeloid.obs['percent.mito'] < 0.15]
ad_myeloid.layers['counts'] = ad_myeloid.X.copy() 
ad_myeloid.layers['raw_counts'] = ad_myeloid.X.copy() 
ad_myeloid.layers['spliced'] = ad_myeloid.X.copy() 
ad_myeloid.layers['unspliced'] = ad_myeloid.X.copy()
ad_myeloid.raw = ad_myeloid

4.2 Preprocessing

We’ll begin our analysis by running our cells through a standard preprocessing pipeline composed of QC, integration, normalization, HVG selection, dimension reduction, & clustering.

4.2.1 QC

We filter out cells with a sequencing depth of less than 1,000, then remove out genes expressed in less than 5 cells.

Code 
sc.pp.filter_cells(ad_myeloid, min_counts=1000)
sc.pp.filter_genes(ad_myeloid, min_cells=5)

4.2.2 HVG selection

Here we select the top 4000 most highly-variable genes (HVGs), which we’ll use as input to PCA, integration, etc.

Code 
sc.pp.highly_variable_genes(
    ad_myeloid, 
    n_top_genes=4000, 
    flavor='seurat_v3', 
    layer='counts',
    subset=True
)

4.2.3 Integration with scVI

Using the scVI library, we perform celltype-aware integration across batches by training a variational autoencoder (VAE) for 150 epochs. This provides us with a 20-dimensional integrated latent space embedding that is (hopefully) free of batch effects.

Code 
scvi.settings.seed = 312
Code 
scvi.settings.num_threads = 16
scvi.model.SCVI.setup_anndata(
    ad_myeloid, 
    layer='counts', 
    batch_key='fetal.ids', 
    labels_key='celltype'
)
int_model = scvi.model.SCVI(
    ad_myeloid, 
    n_layers=2, 
    n_hidden=72, 
    n_latent=20, 
    gene_likelihood='nb', 
    dispersion='gene-label'
)
int_model.train(
    early_stopping=True,
    accelerator='cpu', 
    max_epochs=150
)
Code 
ad_myeloid.obsm['X_scVI'] = int_model.get_latent_representation()
Code 
sc.pl.embedding(
    ad_myeloid, 
    basis='scVI', 
    color='celltype',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5, 
    show=False
)
plt.gca().set_xlabel('scVI 1')
plt.gca().set_ylabel('scVI 2')
plt.show()
Figure 1: scVI embedding colored by celltype.

4.2.4 Normalization

We next depth- and log1p-normalize the mRNA counts matrix.

Code 
ad_myeloid.X = sc.pp.normalize_total(ad_myeloid, target_sum=1e4, inplace=False)['X']
sc.pp.log1p(ad_myeloid)

4.2.5 PCA embedding

Using PCA we generate a 50-dimensional linear reduction of the normalized counts.

Code 
sc.pp.scale(ad_myeloid)
sc.tl.pca(
    ad_myeloid, 
    n_comps=50, 
    random_state=312, 
    use_highly_variable=True
)

Plotting the first two PCs shows clear separation by celltype, though variation in the immature celltypes seems less clear.

Code 
sc.pl.embedding(
    ad_myeloid, 
    basis='pca', 
    color='celltype',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5,
    show=False
)
plt.gca().set_xlabel('PC 1')
plt.gca().set_ylabel('PC 2')
plt.show()
Figure 2: PCA embedding colored by celltype.

4.2.6 SNN graph estimation

Next, we identify the 50 nearest-neighbors (NNs) for each cell using the cosine distance in the scVI latent space.

Code 
sc.pp.neighbors(
    ad_myeloid, 
    n_neighbors=50,
    n_pcs=None,  
    metric='cosine', 
    random_state=312, 
    use_rep='X_scVI'
)

4.2.7 Clustering

Using the Leiden algorithm we partition the SNN graph into clusters.

Code 
sc.tl.leiden(
    ad_myeloid, 
    resolution=0.3, 
    random_state=312
)

The identified clusters roughly match our celltypes.

Code 
sc.pl.embedding(
    ad_myeloid, 
    basis='pca', 
    color='leiden',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5,
    show=False
)
plt.gca().set_xlabel('PC 1')
plt.gca().set_ylabel('PC 2')
plt.show()
Figure 3: PCA embedding colored by Leiden cluster.

4.2.8 Moment-based imputation

Moving on, we create smoothed versions of the counts across each cell’s 30 NNs.

Code 
scv.pp.moments(
    ad_myeloid, 
    n_pcs=None,
    use_rep='X_scVI', 
    n_neighbors=30
)

4.2.9 Undirected PAGA embedding

We use the PAGA algorithm to generate a graph abstraction of the relationships between celltypes.

Code 
sc.tl.paga(ad_myeloid, groups='celltype')

Visualizing the results shows a pretty linear structure that matches what we expect biologically - stem cells at the root, monocytes & macrophages in the middle, and Kupffer cells at the end.

Code 
sc.pl.paga(
    ad_myeloid, 
    fontoutline=True, 
    fontsize=12, 
    frameon=True, 
    random_state=312, 
    show=False
)
plt.gca().set_xlabel('PAGA 1')
plt.gca().set_ylabel('PAGA 2')
plt.show()
Figure 4: PAGA embedding colored by celltype. Lines are weighted by (undirected) connectivity strength.

4.2.10 UMAP embedding

We then create a nonlinear dimension reduction of the data using the UMAP algorithm.

Code 
sc.tl.umap(ad_myeloid, random_state=312)

In UMAP space the transitions between celltypes are clear and well-represented.

Code 
sc.pl.embedding(
    ad_myeloid, 
    basis='umap', 
    color='celltype',
    title='',
    frameon=True, 
    size=30, 
    alpha=0.5,
    show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 5: UMAP embedding colored by celltype.

4.2.11 Force-directed graph embedding

Using the ForceAtlas2 algorithm we generate a 2-dimensional graph layout of our cells.

Code 
sc.tl.draw_graph(
    ad_myeloid, 
    layout='fa', 
    random_state=312
)

The force-directed graph embedding is more visually appealing and cleaner than the UMAP embedding, so we’ll most likely use it instead going forwards.

Code 
sc.pl.draw_graph(
    ad_myeloid, 
    color='celltype', 
    title='', 
    size=30, 
    alpha=0.5, 
    show=False, 
    frameon=True
)
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 6: Force-directed graph embedding colored by celltype.

4.2.12 Diffusion map embedding

Our last embedding will be generated using diffusion maps, which explicitly attempts to preserve transitional structures in the data.

Code 
sc.tl.diffmap(
    ad_myeloid, 
    random_state=312, 
    n_comps=16
)
ad_myeloid.obsm['X_diffmap_old'] = ad_myeloid.obsm['X_diffmap']
ad_myeloid.obsm['X_diffmap'] = ad_myeloid.obsm['X_diffmap'][:, 1:]

While the diffusion map embedding captures some characteristics of the data, it masks variation in the mature macrophage celltypes.

Code 
sc.pl.embedding(
    ad_myeloid, 
    basis='diffmap', 
    color='celltype',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5, 
    components='1,2', 
    show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 7: Diffusion map embedding colored by celltype.

4.2.13 Diffusion pseudotime estimation

We next compute a diffusion pseudotime estimate for each cell using the first diffusion component.

Code 
ad_myeloid.uns['iroot'] = np.argmin(ad_myeloid.obsm['X_diffmap'][:, 0])
sc.tl.dpt(ad_myeloid, n_dcs=1)

Overall the diffusion pseudotime captures myeloid differentiation progression well, but it assigns very similar values to all the mature Kupffer cells which isn’t ideal.

Code 
sc.pl.embedding(
    ad_myeloid, 
    basis='diffmap', 
    color='dpt_pseudotime',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5, 
    components='1,2', 
    show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 8: Diffusion map embedding colored by diffusion pseudotime.

4.3 Save preprocessed data

We save the processed AnnData object to disk.

Code 
ad_myeloid.write_h5ad('../../data/hematopoesis/ad_myeloid.h5ad')

We also save the embeddings, counts matrix, metadata, etc. so that we can create a Seurat object.

Code 
ad_myeloid.obs.reset_index(inplace=False).rename(columns={'index': 'cell'}).to_csv('../../data/hematopoesis/conversion_files/cell_metadata.csv')
ad_myeloid.var.reset_index(inplace=False).rename(columns={'index': 'gene'}).to_csv('../../data/hematopoesis/conversion_files/gene_metadata.csv')
pd.DataFrame(ad_myeloid.obsm['X_draw_graph_fa']).reset_index(inplace=False).to_csv('../../data/hematopoesis/conversion_files/graph_embed.csv')
pd.DataFrame(ad_myeloid.obsm['X_umap']).reset_index(inplace=False).to_csv('../../data/hematopoesis/conversion_files/umap_embed.csv')
pd.DataFrame(ad_myeloid.obsm['X_pca']).reset_index(inplace=False).to_csv('../../data/hematopoesis/conversion_files/pca_embed.csv')
pd.DataFrame(ad_myeloid.obsm['X_scVI']).reset_index(inplace=False).to_csv('../../data/hematopoesis/conversion_files/scvi_embed.csv')
pd.DataFrame(ad_myeloid.obsm['X_diffmap']).reset_index(inplace=False).to_csv('../../data/hematopoesis/conversion_files/diffmap_embed.csv')
mmwrite('../../data/hematopoesis/conversion_files/raw_counts.mtx', a=ad_myeloid.layers['raw_counts'])

5 Analysis

5.1 Read data into R

We start our analysis by importing the preprocessed data from Python into R. Expand the code block for details.

Code 
cell_metadata <- readr::read_csv("../../data/hematopoesis/conversion_files/cell_metadata.csv", 
                                 show_col_types = FALSE, 
                                 col_select = -1) %>% 
                 as.data.frame() %>% 
                 magrittr::set_rownames(.$cell)
gene_metadata <- readr::read_csv("../../data/hematopoesis/conversion_files/gene_metadata.csv", 
                                 show_col_types = FALSE, 
                                 col_select = -1)
embed_PCA <- readr::read_csv("../../data/hematopoesis/conversion_files/pca_embed.csv",
                             show_col_types = FALSE, 
                             col_select = -c(1:2)) %>% 
             as.data.frame() %>% 
             magrittr::set_colnames(paste0("PC", 1:ncol(.))) %>% 
             magrittr::set_rownames(cell_metadata$cell)
embed_PCA <- CreateDimReducObject(embeddings = as.matrix(embed_PCA),
                                  key = "PC_", 
                                  global = TRUE, 
                                  assay = "RNA")
embed_FA <- readr::read_csv("../../data/hematopoesis/conversion_files/graph_embed.csv", 
                            show_col_types = FALSE, 
                            col_select = -c(1:2)) %>% 
            as.data.frame() %>% 
            magrittr::set_colnames(paste0("FA", 1:ncol(.))) %>% 
            magrittr::set_rownames(cell_metadata$cell)
embed_FA <- CreateDimReducObject(embeddings = as.matrix(embed_FA),
                                 key = "FA_", 
                                 global = TRUE, 
                                 assay = "RNA")
embed_UMAP <- readr::read_csv("../../data/hematopoesis/conversion_files/umap_embed.csv", 
                              show_col_types = FALSE, 
                               col_select = -c(1:2)) %>% 
              as.data.frame() %>% 
              magrittr::set_colnames(paste0("UMAP", 1:ncol(.))) %>% 
              magrittr::set_rownames(cell_metadata$cell)
embed_UMAP <- CreateDimReducObject(embeddings = as.matrix(embed_UMAP),
                                   key = "UMAP_", 
                                   global = TRUE, 
                                   assay = "RNA")
embed_diffmap <- readr::read_csv("../../data/hematopoesis/conversion_files/diffmap_embed.csv", 
                                 show_col_types = FALSE, 
                                  col_select = -c(1:2)) %>% 
                 as.data.frame() %>% 
                 magrittr::set_colnames(paste0("DC", 1:ncol(.))) %>% 
                 magrittr::set_rownames(cell_metadata$cell)
embed_diffmap <- CreateDimReducObject(embeddings = as.matrix(embed_diffmap),
                                      key = "DC_", 
                                      global = TRUE, 
                                      assay = "RNA")
embed_scVI <- readr::read_csv("../../data/hematopoesis/conversion_files/scvi_embed.csv", 
                              show_col_types = FALSE, 
                               col_select = -c(1:2)) %>% 
              as.data.frame() %>% 
              magrittr::set_colnames(paste0("scVI", 1:ncol(.))) %>% 
              magrittr::set_rownames(cell_metadata$cell)
embed_scVI <- CreateDimReducObject(embeddings = as.matrix(embed_scVI),
                                   key = "scVI_", 
                                   global = TRUE, 
                                   assay = "RNA")
raw_counts <- Matrix::t(Matrix::readMM("../../data/hematopoesis/conversion_files/raw_counts.mtx"))
colnames(raw_counts) <- cell_metadata$cell
rownames(raw_counts) <- gene_metadata$gene
raw_counts <- CreateAssayObject(counts = raw_counts,
                                min.cells = 0,
                                min.features = 0, 
                                key = "rna")

We can now create a Seurat object, to which we add our various embeddings and metadata. After adding all the various bits and pieces, we normalize the raw counts and re-identify HVGs.

Code 
seu_myeloid <- CreateSeuratObject(raw_counts, 
                                  assay = "RNA",
                                  meta.data = cell_metadata, 
                                  project = "hematopoesis", 
                                  min.cells = 0, 
                                  min.features = 0)
seu_myeloid@meta.data <- mutate(seu_myeloid@meta.data, 
                                celltype = factor(celltype, levels = c("HSC", 
                                                                       "Kupffer cell", 
                                                                       "Monocyte-macrophage", 
                                                                       "Monocyte", 
                                                                       "Monocyte precursor", 
                                                                       "Neutrophil-myeloid progenitor")))
seu_myeloid@reductions$pca <- embed_PCA
seu_myeloid@reductions$diffmap <- embed_diffmap
seu_myeloid@reductions$fa <- embed_FA
seu_myeloid@reductions$umap <- embed_UMAP
seu_myeloid@reductions$scvi <- embed_scVI
seu_myeloid <- NormalizeData(seu_myeloid, verbose = FALSE) %>% 
               FindVariableFeatures(selection.method = "vst", verbose = FALSE)

5.2 Pseudotime estimation

Using Slingshot we identify a single lineage of pseudotime across the force-directed graph embedding.

Code 
sling_res <- slingshot(seu_myeloid@reductions$fa@cell.embeddings[, 1:2], 
                       clusterLabels = as.factor(seu_myeloid$celltype), 
                       start.clus = c("HSC"), 
                       approx_points = 1000)
sling_curves <- slingCurves(sling_res, as.df = TRUE)
sling_mst <- slingMST(sling_res, as.df = TRUE)
sling_pt <- slingPseudotime(sling_res) %>%
            as.data.frame() %>%
            mutate(cell = rownames(.), .before = 1) %>%
            rename(PT = Lineage1) %>% 
            mutate(PT = (PT - min(PT)) / (max(PT) - min(PT)))

Visualizing the results from Slingshot shows us that our biological expectations have been met.

Code 
p0 <- data.frame(Embeddings(seu_myeloid, "fa")) %>% 
      mutate(celltype = seu_myeloid$celltype) %>% 
      ggplot(aes(x = FA_1, y = FA_2, color = celltype)) + 
      geom_point(size = 1.5, 
                 alpha = 0.5, 
                 stroke = 0) + 
      geom_path(data = sling_mst, mapping = aes(x = FA_1, y = FA_2, group = Lineage), 
                linewidth = 1.25, 
                color = "black") + 
      geom_point(data = sling_mst, mapping = aes(x = FA_1, y = FA_2, fill = Cluster), 
                color = "black", 
                shape = 21, 
                size = 4.5, 
                stroke = 1.25, 
                show.legend = FALSE) + 
      scale_color_manual(values = palette_celltype) + 
      scale_fill_manual(values = palette_celltype) + 
      labs(x = "FA 1", y = "FA 2") +
      theme_scLANE(umap = TRUE) +  
      theme(legend.title = element_blank()) + 
      guide_umap()
p1 <- data.frame(Embeddings(seu_myeloid, "fa")) %>% 
      mutate(celltype = seu_myeloid$celltype) %>% 
      ggplot(aes(x = FA_1, y = FA_2, color = celltype)) + 
      geom_point(size = 1.5, 
                 alpha = 0.5, 
                 stroke = 0) + 
      geom_path(data = sling_curves,
                mapping = aes(x = FA_1, y = FA_2, group = Lineage), 
                color = "black", 
                linewidth = 1.5, 
                alpha = 0.5, 
                lineend = "round") + 
      scale_color_manual(values = palette_celltype) + 
      labs(x = "FA 1", y = "FA 2") + 
      theme_scLANE(umap = TRUE) + 
      theme(legend.title = element_blank()) + 
      guide_umap()
p2 <- (p0 / p1) + plot_layout(guides = "collect", axes = "collect")
p2
Figure 9: Force-directed graph embedding overlaid with the MST and principal curves estimated by Slingshot.

5.3 Trajectory DE testing with scLANE

We use the scLANE package to perform trajectory DE testing across pseudotime time for the top 4000 HVGs.

Code 
cell_offset <- createCellOffset(seu_myeloid)
pt_df <- data.frame(PT = sling_pt$PT)
candidate_genes <- chooseCandidateGenes(seu_myeloid, 
                                        group.by.subject = FALSE, 
                                        n.desired.genes = 4000L)
Warning in asMethod(object): sparse->dense coercion: allocating vector of size
1.1 GiB
Code 
scLANE_models <- testDynamic(seu_myeloid, 
                             pt = pt_df, 
                             genes = candidate_genes, 
                             size.factor.offset = cell_offset, 
                             n.cores = 16L,
                             verbose = FALSE)
Warning in asMethod(object): sparse->dense coercion: allocating vector of size
1.1 GiB
scLANE testing in GLM mode completed for 4000 genes across 1 lineage in 2.814 hours
Code 
scLANE_de_res <- getResultsDE(scLANE_models)

5.4 tradeSeq DE analysis

Code 
ts_start <- Sys.time()
bioc_par <- BiocParallel::MulticoreParam(workers = 16L, RNGseed = 312)
RNA_counts <- as.matrix(seu_myeloid@assays$RNA$counts)[candidate_genes, ]
Warning in asMethod(object): sparse->dense coercion: allocating vector of size
1.1 GiB
Code 
k_eval <- evaluateK(RNA_counts, 
                    pseudotime = pt_df, 
                    cellWeights = matrix(rep(1, nrow(pt_df)), ncol = 1), 
                    offset = log(1 / cell_offset), 
                    k = 3:10, 
                    plot = FALSE, 
                    nGenes = 500, 
                    verbose = FALSE, 
                    parallel = TRUE, 
                    BPPARAM = bioc_par)
best_k <- c(3:10)[which.min(abs(colMeans(k_eval - rowMeans(k_eval))))]  # choose k w/ lowest MAD from mean AIC 
ts_models <- fitGAM(RNA_counts, 
                    pseudotime = pt_df, 
                    cellWeights = matrix(rep(1, nrow(pt_df)), ncol = 1), 
                    offset = log(1 / cell_offset), 
                    nknots = best_k, 
                    sce = FALSE, 
                    parallel = TRUE, 
                    verbose = FALSE, 
                    BPPARAM = bioc_par)
BiocParallel::bpstop(bioc_par)
ts_de_res <- associationTest(ts_models, global = TRUE) %>% 
                             arrange(desc(waldStat)) %>% 
                             mutate(gene = rownames(.), 
                                    pvalue_adj = p.adjust(pvalue, method = "fdr"), 
                                    gene_dynamic_overall = if_else(pvalue_adj < 0.01, 1, 0)) %>% 
                             relocate(gene)
ts_end <- Sys.time()
ts_diff <- ts_end - ts_start

The runtime for tradeSeq with 5 knots per gene is:

Code 
ts_diff
Time difference of 2.234959 hours

5.5 Cell fate analysis

Lastly, we’ll perform an analysis of cell fate specification using the CellRank Python package.

5.5.1 CytoTRACE kernel

We begin by computing a CytoTRACE kernel; CytoTRACE uses the number of genes expressed in a cell as a proxy measurement for differentiation potential, with the assumption being that cells expressing more genes are more likely to be progenitor celltypes & vice versa.

Code 
ctk = CytoTRACEKernel(ad_myeloid).compute_cytotrace()

We plot the CytoTRACE scores below. Interestingly, in addition to the expected result of HSCs having high CytoTRACE scores we see a subcluster of Kupffer cells with high scores as well.

Code 
sc.pl.embedding(
    ad_myeloid,
    basis='draw_graph_fa',
    color=['ct_score', 'ct_pseudotime'],
    color_map='magma', 
    show=False, 
    size=30, 
    alpha=0.5, 
    title=['CytoTRACE score', 'CytoTRACE Pseudotime']
)
plt.show()
Figure 10: Force-directed graph embedding colored by CytoTRACE score and pseudotime. Higher scores indicate higher differentiation potential.

Next, we estimate a cell-cell transition probability matrix.

Code 
ctk.compute_transition_matrix(threshold_scheme='soft', show_progress_bar=False)

Plotting the projection of that matrix on our force-directed graph embedding shows messy / reversed directionality. This indicates that CytoTRACE might not work very well on this dataset, and at the very least should be augmented with another data source.

Code 
ctk.plot_projection(
  basis='draw_graph_fa', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.5, 
  linewidth=2, 
  frameon=True, 
  show=False
)
/blue/rbacher/j.leary/py_envs/scLANE_env2/lib/python3.10/site-packages/scvelo/plotting/utils.py:63: DeprecationWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, pd.CategoricalDtype) instead
  return isinstance(c, str) and c in data.obs.keys() and cat(data.obs[c])
/blue/rbacher/j.leary/py_envs/scLANE_env2/lib/python3.10/site-packages/scvelo/plotting/utils.py:63: DeprecationWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, pd.CategoricalDtype) instead
  return isinstance(c, str) and c in data.obs.keys() and cat(data.obs[c])
Code 
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 11: Force-directed graph embedding streamline plot showing differentiation directions as inferred from CytoTRACE scores.

5.5.2 Pseudotime kernel

Next, we use the pseudotime estimates from Slingshot to create a pseudotime-based kernel and estimate another cell-cell transition probability matrix.

Code 
ad_myeloid.obs['sling_PT'] = r.sling_pt['PT'].tolist()
pk = PseudotimeKernel(ad_myeloid, time_key='sling_PT')
pk.compute_transition_matrix(threshold_scheme='soft', show_progress_bar=False)

Plotting the projection of the pseudotime kernel’s transition probability matrix shows a much more reasonable portrayal of differentiation directionality.

Code 
pk.plot_projection(
  basis='draw_graph_fa', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.75, 
  linewidth=2, 
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 12: Force-directed graph embedding overlaid with streamline arrows derived from Slingshot pseudotime.

5.5.3 Combined kernel

We combine the two kernel using unequal weighting, with the pseudotime kernel being given priority due to its better overall capture of differentiation.

Code 
ck = 0.2 * ctk + 0.8 * pk

We visualize the combined kernel’s projection below. It seems to reliably portray the direction of differentiation in our dataset.

Code 
ck.plot_projection(
  basis='draw_graph_fa', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.75, 
  linewidth=2, 
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 13: Force-directed graph embedding overlaid with streamline arrows derived from the combined kernel.

We can also simulate random walks on the high-dimensional cell-cell graph starting from the HSC population. We see that nearly all of the random walks terminate in the Kupffer cell population.

Code 
ck.plot_random_walks(
    n_sims=200,
    start_ixs={'celltype': 'HSC'},
    basis='draw_graph_fa',
    color='celltype',
    legend_loc='right margin',
    seed=312, 
    frameon=True, 
    title='', 
    linewidth=0.5, 
    linealpha=0.25, 
    size=30, 
    alpha=0.5, 
    ixs_legend_loc='lower', 
    show_progress_bar=False
)
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 14: Force-directed graph embedding overlaid with random walks along the cell-cell graph derived from the combined kernel. Starting points are highlighted in black and ending points in yellow.

6 Save data

Finally, we save our Seurat object, the output from Slingshot, and the models from scLANE and tradeSeq.

Code 
saveRDS(seu_myeloid, file = "../../data/hematopoesis/seu_myeloid.Rds")
saveRDS(sling_res, file = "../../data/hematopoesis/sling_res.Rds")
saveRDS(scLANE_models, file = "../../data/hematopoesis/scLANE_models.Rds")
saveRDS(ts_models, file = "../../data/hematopoesis/ts_models.Rds")

7 Session info

Code 
sessioninfo::session_info()
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.3 (2024-02-29)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/New_York
 date     2024-10-11
 pandoc   2.9.2.1 @ /usr/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package              * version    date (UTC) lib source
 abind                  1.4-8      2024-09-12 [2] CRAN (R 4.3.3)
 AnnotationDbi          1.62.2     2023-07-02 [2] Bioconductor
 backports              1.5.0      2024-05-23 [1] CRAN (R 4.2.3)
 bigassertr             0.1.6      2023-01-10 [1] CRAN (R 4.3.2)
 bigparallelr           0.3.2      2021-10-02 [1] CRAN (R 4.3.2)
 bigstatsr              1.5.12     2022-10-14 [1] CRAN (R 4.3.2)
 Biobase              * 2.60.0     2023-04-25 [2] Bioconductor
 BiocFileCache          2.13.1     2024-10-05 [1] Github (Bioconductor/BiocFileCache@30abbe4)
 BiocGenerics         * 0.46.0     2023-04-25 [2] Bioconductor
 BiocParallel           1.34.2     2023-05-22 [2] Bioconductor
 biomaRt              * 2.56.1     2023-06-09 [1] Bioconductor
 Biostrings             2.68.1     2023-05-16 [2] Bioconductor
 bit                    4.5.0      2024-09-20 [2] CRAN (R 4.3.3)
 bit64                  4.5.2      2024-09-22 [2] CRAN (R 4.3.3)
 bitops                 1.0-8      2024-07-29 [2] CRAN (R 4.3.3)
 blob                   1.2.4      2023-03-17 [2] CRAN (R 4.3.1)
 boot                   1.3-31     2024-08-28 [4] CRAN (R 4.3.3)
 broom                  1.0.6      2024-05-17 [1] CRAN (R 4.2.3)
 broom.mixed            0.2.9.5    2024-04-01 [1] CRAN (R 4.2.3)
 cachem                 1.1.0      2024-05-16 [2] CRAN (R 4.3.3)
 cli                    3.6.2      2023-12-11 [1] CRAN (R 4.2.3)
 cluster                2.1.6      2023-12-01 [4] CRAN (R 4.3.2)
 coda                   0.19-4.1   2024-01-31 [2] CRAN (R 4.3.2)
 codetools              0.2-20     2024-03-31 [4] CRAN (R 4.3.3)
 colorspace             2.1-1      2024-07-26 [2] CRAN (R 4.3.3)
 cowplot                1.1.3      2024-01-22 [2] CRAN (R 4.3.3)
 crayon                 1.5.3      2024-06-20 [2] CRAN (R 4.3.3)
 curl                   5.2.3      2024-09-20 [2] CRAN (R 4.3.3)
 data.table             1.16.0     2024-08-27 [2] CRAN (R 4.3.3)
 DBI                    1.2.3      2024-06-02 [2] CRAN (R 4.3.3)
 dbplyr                 2.5.0      2024-03-19 [1] CRAN (R 4.3.3)
 DelayedArray           0.26.7     2023-07-28 [2] Bioconductor
 DelayedMatrixStats     1.22.6     2023-08-28 [2] Bioconductor
 deldir                 2.0-4      2024-02-28 [2] CRAN (R 4.3.3)
 digest                 0.6.35     2024-03-11 [1] CRAN (R 4.2.3)
 doParallel             1.0.17     2022-02-07 [1] CRAN (R 4.2.3)
 doSNOW                 1.0.20     2022-02-04 [1] CRAN (R 4.3.2)
 dotCall64              1.1-1      2023-11-28 [2] CRAN (R 4.3.1)
 dplyr                * 1.1.4      2023-11-17 [1] CRAN (R 4.2.3)
 edgeR                  3.42.4     2023-05-31 [2] Bioconductor
 emmeans                1.10.2     2024-05-20 [1] CRAN (R 4.3.3)
 estimability           1.5.1      2024-05-12 [1] CRAN (R 4.3.3)
 evaluate               0.24.0     2024-06-10 [1] CRAN (R 4.3.3)
 fansi                  1.0.6      2023-12-08 [2] CRAN (R 4.3.1)
 farver                 2.1.2      2024-05-13 [1] CRAN (R 4.2.3)
 fastDummies            1.7.4      2024-08-16 [2] CRAN (R 4.3.3)
 fastmap                1.2.0      2024-05-15 [2] CRAN (R 4.3.3)
 ff                     4.0.12     2024-01-12 [1] CRAN (R 4.3.2)
 filelock               1.0.3      2023-12-11 [2] CRAN (R 4.3.3)
 fitdistrplus           1.2-1      2024-07-12 [2] CRAN (R 4.3.3)
 flock                  0.7        2016-11-12 [1] CRAN (R 4.3.2)
 forcats                1.0.0      2023-01-29 [1] CRAN (R 4.3.2)
 foreach                1.5.2      2022-02-02 [1] CRAN (R 4.2.3)
 furrr                  0.3.1      2022-08-15 [1] CRAN (R 4.3.2)
 future                 1.33.2     2024-03-26 [1] CRAN (R 4.2.3)
 future.apply           1.11.2     2024-03-28 [2] CRAN (R 4.3.3)
 gamlss                 5.4-22     2024-03-20 [1] CRAN (R 4.3.2)
 gamlss.data            6.0-6      2024-03-14 [1] CRAN (R 4.3.2)
 gamlss.dist            6.1-1      2023-08-23 [1] CRAN (R 4.3.2)
 geeM                   0.10.1     2018-06-18 [1] CRAN (R 4.3.2)
 generics               0.1.3      2022-07-05 [1] CRAN (R 4.2.0)
 GenomeInfoDb         * 1.36.4     2023-10-02 [2] Bioconductor
 GenomeInfoDbData       1.2.10     2023-07-26 [2] Bioconductor
 GenomicRanges        * 1.52.1     2023-10-08 [2] Bioconductor
 ggplot2              * 3.5.1      2024-04-23 [1] CRAN (R 4.2.3)
 ggrepel                0.9.2      2022-11-06 [1] CRAN (R 4.2.2)
 ggridges               0.5.6      2024-01-23 [2] CRAN (R 4.3.3)
 glm2                 * 1.2.1      2018-08-11 [1] CRAN (R 4.3.2)
 glmmTMB                1.1.9      2024-03-20 [1] CRAN (R 4.3.3)
 globals                0.16.3     2024-03-08 [2] CRAN (R 4.3.3)
 glue                   1.8.0      2024-09-30 [2] CRAN (R 4.3.3)
 goftest                1.2-3      2021-10-07 [2] CRAN (R 4.3.1)
 gridExtra              2.3        2017-09-09 [2] CRAN (R 4.3.1)
 gtable                 0.3.5      2024-04-22 [2] CRAN (R 4.3.3)
 hms                    1.1.3      2023-03-21 [2] CRAN (R 4.3.1)
 htmltools              0.5.8.1    2024-04-04 [1] CRAN (R 4.2.3)
 htmlwidgets            1.6.4      2023-12-06 [2] CRAN (R 4.3.3)
 httpuv                 1.6.15     2024-03-26 [2] CRAN (R 4.3.3)
 httr                   1.4.7      2023-08-15 [2] CRAN (R 4.3.3)
 ica                    1.0-3      2022-07-08 [2] CRAN (R 4.3.1)
 igraph                 2.0.3      2024-03-13 [1] CRAN (R 4.2.3)
 IRanges              * 2.34.1     2023-06-22 [2] Bioconductor
 irlba                  2.3.5.1    2022-10-03 [2] CRAN (R 4.3.1)
 iterators              1.0.14     2022-02-05 [2] CRAN (R 4.3.1)
 jsonlite               1.8.8      2023-12-04 [1] CRAN (R 4.2.3)
 KEGGREST               1.40.1     2023-09-29 [2] Bioconductor
 KernSmooth             2.23-24    2024-05-17 [4] CRAN (R 4.3.3)
 knitr                  1.48       2024-07-07 [2] CRAN (R 4.3.3)
 labeling               0.4.3      2023-08-29 [2] CRAN (R 4.3.1)
 later                  1.3.2      2023-12-06 [2] CRAN (R 4.3.3)
 lattice                0.22-6     2024-03-20 [4] CRAN (R 4.3.3)
 lazyeval               0.2.2      2019-03-15 [2] CRAN (R 4.3.1)
 leiden                 0.4.3.1    2023-11-17 [2] CRAN (R 4.3.1)
 lifecycle              1.0.4      2023-11-07 [2] CRAN (R 4.3.1)
 limma                  3.56.2     2023-06-04 [2] Bioconductor
 listenv                0.9.1      2024-01-29 [2] CRAN (R 4.3.3)
 lme4                   1.1-35.5   2024-07-03 [2] CRAN (R 4.3.3)
 lmtest                 0.9-40     2022-03-21 [2] CRAN (R 4.3.1)
 locfit                 1.5-9.10   2024-06-24 [2] CRAN (R 4.3.3)
 magrittr             * 2.0.3      2022-03-30 [1] CRAN (R 4.2.0)
 MASS                   7.3-60     2023-05-04 [4] CRAN (R 4.3.1)
 Matrix                 1.6-5      2024-01-11 [1] CRAN (R 4.2.3)
 MatrixGenerics       * 1.12.3     2023-07-30 [2] Bioconductor
 matrixStats          * 1.4.1      2024-09-08 [2] CRAN (R 4.3.3)
 memoise                2.0.1      2021-11-26 [2] CRAN (R 4.3.1)
 mgcv                   1.9-1      2023-12-21 [4] CRAN (R 4.3.2)
 mime                   0.12       2021-09-28 [2] CRAN (R 4.3.1)
 miniUI                 0.1.1.1    2018-05-18 [2] CRAN (R 4.3.1)
 minqa                  1.2.7      2024-05-20 [1] CRAN (R 4.2.3)
 munsell                0.5.1      2024-04-01 [2] CRAN (R 4.3.3)
 mvtnorm                1.3-1      2024-09-03 [2] CRAN (R 4.3.3)
 nlme                   3.1-166    2024-08-14 [4] CRAN (R 4.3.3)
 nloptr                 2.1.1      2024-06-25 [2] CRAN (R 4.3.3)
 numDeriv               2016.8-1.1 2019-06-06 [2] CRAN (R 4.3.1)
 paletteer              1.6.0      2024-01-21 [1] CRAN (R 4.2.3)
 parallelly             1.37.1     2024-02-29 [1] CRAN (R 4.2.3)
 patchwork            * 1.2.0      2024-01-08 [1] CRAN (R 4.2.3)
 pbapply                1.7-2      2023-06-27 [2] CRAN (R 4.3.1)
 pillar                 1.9.0      2023-03-22 [2] CRAN (R 4.3.1)
 pkgconfig              2.0.3      2019-09-22 [2] CRAN (R 4.3.1)
 plotly                 4.10.4     2024-01-13 [2] CRAN (R 4.3.3)
 plyr                   1.8.9      2023-10-02 [2] CRAN (R 4.3.3)
 png                    0.1-8      2022-11-29 [2] CRAN (R 4.3.1)
 polyclip               1.10-7     2024-07-23 [2] CRAN (R 4.3.3)
 prettyunits            1.2.0      2023-09-24 [2] CRAN (R 4.3.3)
 princurve            * 2.1.6      2021-01-18 [1] CRAN (R 4.3.2)
 prismatic              1.1.2      2024-04-10 [1] CRAN (R 4.2.3)
 progress               1.2.3      2023-12-06 [2] CRAN (R 4.3.3)
 progressr              0.14.0     2023-08-10 [2] CRAN (R 4.3.2)
 promises               1.3.0      2024-04-05 [2] CRAN (R 4.3.3)
 ps                     1.8.0      2024-09-12 [2] CRAN (R 4.3.3)
 purrr                  1.0.2      2023-08-10 [1] CRAN (R 4.2.3)
 R6                     2.5.1      2021-08-19 [2] CRAN (R 4.3.1)
 RANN                   2.6.2      2024-08-25 [2] CRAN (R 4.3.3)
 rappdirs               0.3.3      2021-01-31 [2] CRAN (R 4.3.1)
 RColorBrewer           1.1-3      2022-04-03 [2] CRAN (R 4.3.1)
 Rcpp                   1.0.13     2024-07-17 [2] CRAN (R 4.3.3)
 RcppAnnoy              0.0.22     2024-01-23 [2] CRAN (R 4.3.3)
 RcppEigen              0.3.4.0.0  2024-02-28 [1] CRAN (R 4.2.3)
 RcppHNSW               0.6.0      2024-02-04 [2] CRAN (R 4.3.3)
 RCurl                  1.98-1.16  2024-07-11 [2] CRAN (R 4.3.3)
 readr                  2.1.5      2024-01-10 [2] CRAN (R 4.3.3)
 rematch2               2.1.2      2020-05-01 [2] CRAN (R 4.3.1)
 reshape2               1.4.4      2020-04-09 [2] CRAN (R 4.3.1)
 reticulate           * 1.35.0     2024-01-31 [1] CRAN (R 4.2.3)
 rlang                  1.1.4      2024-06-04 [1] CRAN (R 4.3.3)
 rmarkdown              2.27       2024-05-17 [1] CRAN (R 4.3.3)
 rmio                   0.4.0      2022-02-17 [1] CRAN (R 4.3.2)
 ROCR                   1.0-11     2020-05-02 [2] CRAN (R 4.3.1)
 RSpectra               0.16-2     2024-07-18 [2] CRAN (R 4.3.3)
 RSQLite                2.3.7      2024-05-27 [2] CRAN (R 4.3.3)
 Rtsne                  0.17       2023-12-07 [2] CRAN (R 4.3.3)
 S4Arrays               1.2.1      2024-03-04 [1] Bioconductor 3.18 (R 4.3.2)
 S4Vectors            * 0.38.2     2023-09-22 [2] Bioconductor
 scales                 1.3.0      2023-11-28 [1] CRAN (R 4.2.3)
 scattermore            1.2        2023-06-12 [2] CRAN (R 4.3.1)
 scLANE               * 0.8.4      2024-10-10 [1] Github (jr-leary7/scLANE@bc239c0)
 sctransform            0.4.1      2023-10-19 [1] CRAN (R 4.2.3)
 sessioninfo            1.2.2      2021-12-06 [2] CRAN (R 4.3.1)
 Seurat               * 5.1.0      2024-05-10 [1] CRAN (R 4.3.3)
 SeuratObject         * 5.0.2      2024-05-08 [1] CRAN (R 4.3.3)
 shiny                  1.9.1      2024-08-01 [2] CRAN (R 4.3.3)
 SingleCellExperiment * 1.22.0     2023-04-25 [2] Bioconductor
 slingshot            * 2.6.0      2022-11-01 [1] Bioconductor
 snow                   0.4-4      2021-10-27 [2] CRAN (R 4.3.1)
 sp                   * 2.1-4      2024-04-30 [3] CRAN (R 4.3.3)
 spam                   2.10-0     2023-10-23 [2] CRAN (R 4.3.1)
 sparseMatrixStats      1.12.2     2023-07-02 [2] Bioconductor
 spatstat.data          3.1-2      2024-06-21 [2] CRAN (R 4.3.3)
 spatstat.explore       3.3-2      2024-08-21 [2] CRAN (R 4.3.3)
 spatstat.geom          3.3-3      2024-09-18 [2] CRAN (R 4.3.3)
 spatstat.random        3.3-2      2024-09-18 [2] CRAN (R 4.3.3)
 spatstat.sparse        3.1-0      2024-06-21 [2] CRAN (R 4.3.3)
 spatstat.univar        3.0-1      2024-09-05 [2] CRAN (R 4.3.3)
 spatstat.utils         3.1-0      2024-08-17 [2] CRAN (R 4.3.3)
 stringi                1.8.4      2024-05-06 [2] CRAN (R 4.3.3)
 stringr                1.5.1      2023-11-14 [2] CRAN (R 4.3.2)
 SummarizedExperiment * 1.30.2     2023-06-06 [2] Bioconductor
 survival               3.7-0      2024-06-05 [2] CRAN (R 4.3.3)
 tensor                 1.5        2012-05-05 [2] CRAN (R 4.3.1)
 tibble                 3.2.1      2023-03-20 [2] CRAN (R 4.3.1)
 tidyr                  1.3.1      2024-01-24 [1] CRAN (R 4.2.3)
 tidyselect             1.2.1      2024-03-11 [2] CRAN (R 4.3.3)
 TMB                    1.9.11     2024-04-03 [1] CRAN (R 4.3.3)
 tradeSeq             * 1.14.0     2023-04-25 [1] Bioconductor
 TrajectoryUtils      * 1.6.0      2022-11-01 [1] Bioconductor
 tzdb                   0.4.0      2023-05-12 [2] CRAN (R 4.3.1)
 utf8                   1.2.4      2023-10-22 [2] CRAN (R 4.3.1)
 uwot                   0.2.2      2024-04-21 [2] CRAN (R 4.3.3)
 vctrs                  0.6.5      2023-12-01 [2] CRAN (R 4.3.1)
 viridis                0.6.5      2024-01-29 [2] CRAN (R 4.3.3)
 viridisLite            0.4.2      2023-05-02 [2] CRAN (R 4.3.1)
 vroom                  1.6.5      2023-12-05 [2] CRAN (R 4.3.3)
 withr                  3.0.0      2024-01-16 [1] CRAN (R 4.2.3)
 xfun                   0.44       2024-05-15 [1] CRAN (R 4.3.3)
 XML                    3.99-0.17  2024-06-25 [2] CRAN (R 4.3.3)
 xml2                   1.3.6      2023-12-04 [2] CRAN (R 4.3.3)
 xtable                 1.8-4      2019-04-21 [2] CRAN (R 4.3.1)
 XVector                0.40.0     2023-04-25 [2] Bioconductor
 yaml                   2.3.10     2024-07-26 [2] CRAN (R 4.3.3)
 zlibbioc               1.46.0     2023-04-25 [2] Bioconductor
 zoo                    1.8-12     2023-04-13 [2] CRAN (R 4.3.1)

 [1] /home/j.leary/r_packages_default
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

─ Python configuration ───────────────────────────────────────────────────────
 python:         /blue/rbacher/j.leary/py_envs/scLANE_env2/bin/python
 libpython:      /blue/rbacher/j.leary/py_envs/scLANE_env2/lib/libpython3.10.so
 pythonhome:     /blue/rbacher/j.leary/py_envs/scLANE_env2:/blue/rbacher/j.leary/py_envs/scLANE_env2
 version:        3.10.14 | packaged by conda-forge | (main, Mar 20 2024, 12:45:18) [GCC 12.3.0]
 numpy:          /home/j.leary/.local/lib/python3.10/site-packages/numpy
 numpy_version:  1.24.2
 
 NOTE: Python version was forced by use_python() function

──────────────────────────────────────────────────────────────────────────────